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Abstract. We present a systematic study of moment evolution in multidimensional stochastic difference 
systems, focusing on characterizing systems whose low-order moments diverge in the neighborhood of a 
stable fixed point. We consider systems with a simple, dominant eigenvalue and stationary, white noise. 
When the noise is small, we obtain general expressions for the approximate asymptotic distribution and 
moment Lyapunov exponents. In the case of larger noise, the second moment is calculated using a different 
approach, which gives an exact result for some types of noise. We analyze the dependence of the moments on 
the system's dimension, relevant system properties, the form of the noise, and the magnitude of the noise. 
We determine a critical value for noise strength, as a function of the unperturbed system's convergence 
rate, above which the second moment diverges and large fluctuations are likely. Analytical results are 
validated by numerical simulations. We show that our results cannot be extended to the continuous time 
limit except in certain special cases. 

PACS. 02.50.Ey Stochastic processes - 02.50.Sk Multivariate analysis - 05.45.Ca Noise 
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i ■ 1 Introduction 
03 1.1 Motivation and previous work 
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The stability of fixed points in a multidimensional system 
is easily ascertained when the system is perfectly deter- 
ministic by using linear stability analysispp. Many real- 
world systems, however, are not perfectly deterministic 
because their interactions are subject to noise |2J. It is 
therefore of interest to consider the effect of a multiplica- 
tive noise term on a linearized system: 



x* = (A + B 



t\ x t- 



(1) 



In this paper we analyze the effect of white, station- 
ary mean noise in discrete systems. This type of noise 
has no effect on a system's stability in mean, because the 
expected value evolves exactly as if the system were un- 
perturbed f ii2.1|) . However, multiplicative noise processes 
cause fluctuations which can be large even if the fixed 
point is stable (figure Q), knocking the system out of the 
linear regime and coupling it to nonlinearities. Even for 
exact linear models, large fluctuations can cause long de- 
lays in convergence. An example of fluctuations in such a 
system is shown is figure ^ where the dotted line shows 
evolution of the first component of x without noise, and 
the solid line shows one instance of evolution with noise. 

Fluctuations in a stochastic system are studied by way 
of the system's moments The pth moment of a multi- 
variate system is simply the expected value of |x| p ; large 
moments, especially the low order moments such as the 
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Fig. 1. Example of fluctuations in a linear system 



second and third, indicate that a system attains large val- 
ues with non- negligible probability For example, the 
system of figure^ has a divergent moments for p > 3. 

Multiplicative noise causes fluctuations because its ef- 
fect is to cause the moments of a system to diverge, even 
when the system converges in mean 315. . In particular, di- 
vergent low-order moments in the neighborhood of a sta- 
ble fixed point are likely to cause the large fluctuations 
described above. The evolution of the moments is thus an 
important consideration in regards to fixed point stability 
in systems whose interactions are subject to noise. 
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The asymptotic behavior of a random system and its 
moments is characterized by the system's Lyapunov expo- 
nent and moment Lyapunov exponents 0. Calculation of 
Lyapunov exponents for multivariate systems is very diffi- 
cult in general, even in simple cases [7J- Stability analysis 
and calculation of Lyapunov exponents for discrete linear 
stochastic systems and random matrix products has been 
a major area of research in mathematics |81l9lllUlllllll21 
m^TH] . control theory [Tall 1 611 17lHH) . physics [T§II2U| . engi- 
neermg mechanics |!H1l22ll23| . and biology 0E1, among 
others. However, the subject of most research has been 
stability in mean, not stability of the moments. The tra- 
ditional approach to determining convergence in random 
systems is to use bounds (above mathematics references; 
see also |25U2fij . for example, for continuous systems). 

There has been little previous work on calculating ex- 
act expressions for moment evolution in multidimensional 
stochastic systems. One exception is [22] for continuous 
2-dimensional systems, the results of which are discussed 
in context in §3 

1.2 Problem statement and notation 

We are studying a system evolving according to the dif- 
ference equation 

x f = (A + B^x*- 1 , (2) 

or 

t 

x^fQtA + imx . (3) 

T = l 

Here x is the system state, a vector of random variables 
and B* is a matrix of white noise processes with mean 0. 
(That is, (Bjj) = and (B^Bfj) ~ 6 tt >.) The initial state 
x° of the system is assumed to be fixed. The eigenvalues 
of the matrix A are A^; the largest eigenvalue Ai or simply 
A is simple 1 and dominant, that is, A > Ai for all i ^ 1. 

The system size is n. We define the mean of the 
to be a, and the variance to be o~\. In the mean value 
approximation, A ss aG where G is the matrix whose ele- 
ments are all 1. We will be diagonalizing A into the form 
PAP -1 , where A is the diagonal matrix of eigenvalues of 
A, and 

u = Pji = normalized right eigenvector (4) 

corresponding to A; 
v = (P^ 1 ) 7 " = non-normalized left eigenvector (5) 
corresponding to A. 

so that |u| = 1. Note that v ■ u = 1 and so |v| > I. 

A vector x* converges in mean if (cc*) converges for 
all i . We use any typical definition for convergence. The 
system's fixed point is stable whenever the system con- 
verges in mean, because the initial state is irrelevant to 
convergence [2B1- 

1 Simple eigenvalues have algebraic multiplicity 1 and thus 
only one associated eigenvector. 



We define the pth moment of the system to be (|x*| p ). 
Moment convergence can be elegantly expressed in terms 
of moment Lyapunov exponents, discussed in and de- 
fined as 



The asymptotic behavior of the pth moment is then given 
by (|x*| p ) ~ e tLp and the pth moment converges if 

L p < 0. 

Finally, in the case that all the elements of the noise 
matrix have the same variance b 2 , we define the critical 
value b 2 to be the level of noise above which the second 
moment diverges. 

1.3 Overview of results 

The central results of this paper are the approximations 
and exact expressions describing the evolution of moments 
of the system (j2J, in particular the second moment. The 
small noise case is treated first using a perturbation ap- 
proach. This approach allows us to calculate the system's 
approximate asymptotic distribution and moment Lya- 
punov exponents. For larger noises, an iteration technique 
is presented which gives both small and large noise results 
for the second moment. For certain types of noise, the it- 
eration method allows us to calculate the second moment 
Lyapunov exponent exactly in any system. These results 
appear to be the first general analytic results for the Lya- 
punov exponents of discrete multivariate systems. 

The analysis of this paper is valid in discrete systems 
with a simple, dominant eigenvalue. The eigenvalue re- 
quirement is satisfied by all nonnegative systems (see ap- 
pendix [bJ) and many arbitrary systems. Nonnegative |29| 
and positive P2KI1 discrete systems arise in Markov mod- 
els, and the fields of biology, population models, economics 
(input-output models), finance, and cooperative problem 
solving, among others. Applications to arbitrary systems 
are too numerous to list. 

Particular results of this paper are as follows. First, we 
show that in the small noise regime, the problem of ap- 
proximating the asymptotic probability distribution of a 
multidimensional system reduces to the scalar case, which 
is trivial ( 12.21 ij4.ll H4.3|l . We thus obtain the expression 

(|xr) = |(x t )| p e te2£i2 ? ii+0 ( £4 ) (7) 

where (x*) is the expected (unperturbed) value of the sys- 
tem at time t, and e is a small parameter which depends 
on the noise and is calculated for various forms of noise 
(equation (|29|1 and table 0. This approximation is justi- 
fied by simulation (figures EdHJl and its accuracy is dis- 
cussed briefly in ^4.51 

In the case of larger noise, the iteration approach of 
section ^presents a methodology for calculating the sec- 
ond moment Lyapunov exponent to any degree of accu- 
racy in any system, provided the noise elements have the 
same variance. The exact value of the Lyapunov exponent 
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is expressed as the largest eigenvalue of a matrix and its 
accuracy is justified in the simulation of figure 

It is shown that the results of the iteration technique 
agree with JJJ for small noise, and with the A — > limit 
for large noise. It is also shown that all results agree with 
the trivial scalar case discussed in H2.2I 

While the unperturbed value of the system depends 
only on the initial state and the dominant eigenvalue in 
the asymptotic limit, the moments depends on other prop- 
erties of the system including the system size and the form 
of the noise. It is shown that 

— the effect of a given level of noise can be magnified, in 
some cases greatly, if the dominant eigenvalue of the 
unperturbed system is ill-conditioned (^2 M.3fl : 

— the destabilizing effect of the noise is damped as the 
number of independent components of noise increases 
( "destructive interference" of independent noises) ( MA\ 
figure 0) ; 

— the destructive interference of independent noises is 
maximized in the mean value limit f HI.2|) and is mit- 
igated by any deviation from this limit f H4.4l figure 

o§. 

— large noise ( H6.2|l . or small noise in systems with a 
very ill-conditioned dominant eigenvalue ( ^5.21 figure 
EJ, almost certainly destabilizes the system. 

We also present a discussion of the critical value for the 
noise variance above which the second moment diverges 
and fluctuations become a major consideration. This dis- 
cussion is largely restricted to the case in which all the 
noise elements have the same variance for simplicity. We 
obtain the following expression for the critical value: 

6 <= = „fc i ~E> Za y ' ^ 

where f u and f v are parameters related to A and the 
type of noise considered which are very close to I in the 
large majority of systems, and k equals I if all the noise 
elements are independent and 2 if they are all correlated. 
This expression is shown to be accurate in both the small 
and large noise cases, and is used to create a stability 
diagram for the system in figures ITT1 and IT51 

The dependence of the critical value on system param- 
eters is discussed. We show that 

— for small noise, the critical value depends weakly on the 
system size and type of noise considered fi Jtj.41 figure 

— for large noise, the critical value depends strongly on 
the system size and type of noise f i|6.4[) : 

— the critical value provides a much more accurate in- 
dication of the level of noise below which the second 
moment converges than a simple bound on convergence 
(123 appendix 0; 

— for most convergent systems subject to small noise, 
the low-order moments diverge only if the unperturbed 
system converges slowly f iJ6.1|) . 

This last statement is especially true for positive systems 
(figure ; note that systems with slow convergence may 



have other problems besides fluctuations due to noise, such 
as large transient behavior [d'2\. 

Finally, we consider the continuous limit and show that 
our results only extend to this limit in certain very special 
cases (SO 

1.4 Paper organization 

The paper is organized as follows. In S|2]we present simple 
preliminary results: asymptotic expressions for the sys- 
tem's expected value, moments in the scalar (n = 1) 
case, and second moment in the case that A = 0. Section 
01 discusses the important properties of the multivariate 
system. Moment evolution is calculated for multivariate 
systems in the small noise limit using a perturbation ap- 
proach in an d the result is discussed. Section [S] uses an 
iteration approach to treat the second moment's evolu- 
tion in the case of larger noise. The critical value of noise 
for second moment divergence is the subject of ^ The 
accuracy of the approximations is justified in numerical 
simulations throughout the paper. Finally, jJ3 presents a 
discussion of the continuous time limit. 



2 Preliminaries 

This section presents a calculation of the expected value 
of the system, as well as calculations of two limiting cases: 
a scalar stochastic system, and a multivariate system with 
no fixed part (noise only). 

A discussion of the expected value of the system and 
its convergence properties is a necessary preliminary step 
to any study of the moments. The computation is trivial, 
as we show, because the noise is white with mean zero. 

Calculation of the moments of a scalar system provides 
a framework which we will apply in the small noise limit of 
the multivariate case (SQJ. The scalar system also provides 
a demonstration of how multiplicative noise leads to a log- 
normal distribution and moment divergence. 

The noise only, A = multivariate system is a system 
in which the moments can be found exactly, yielding the 
zeroth order term for the large- noise limit (see H5.3|l . The 
calculations involved also provide a useful preview of those 
in El 

All of the calculations in this section are quite straight- 
forward and have very likely been presented, in whole or 
part, in some previous work. However, we did not find a 
specific reference with the exception of |18| for a cursory 
treatment of the scalar case. 



2.1 Expected value and unperturbed system 

The expected (average) state of the system and the state 
of the unperturbed system are equivalent since the noise 
is white with mean 0. White noise means that x* _1 and 
B* are independent, so that 

(x t )-(A+(B t ))(x t - 1 )-A(x t - 1 ), 
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since the mean is 0. Thus 



= A*x° = : 



^unperturbed' 



In systems with a simple dominant eigenvalue A, the asymp- 
totic behavior of unperturbed system is completely deter- 
mined by the largest eigenvalue of A 28 . For large t, 



x 

unperturbed 



A'(vx°)u 



(9) 



and the moment Lyapunov exponents are simply 

L° p =p\og\. (10) 

The system will converge to for any initial conditions if 
A < 1, and it will diverge if A i > 1. In the case of stochastic 
matrices with A = 1 , the above formula is accurate because 
A is simple. We are not interested in the case in which x° 
is orthogonal to v. 



2.2 Scalar stochastic system 

In the case n = 1 it is not difficult to determine the asymp- 
totic distribution, as well as the exact expressions for any 
moment. We go through a derivation here because this 
analysis will apply to the small noise multivariate case. 
The scalar system is 



Xt = x J|(a + M> 



where (b T ) — and (b T b T i) = b 2 S rT /. Notice that we ex- 
press time as a subscript in this section, whereas in the 
multidimensional treatment time is a superscript. 



2.2.1 Exact expressions 

We have 

(x p t ) = x p ([a+b T ry 

~ " x ' . (ii) 

In particular, (x t ) — a* and (x 2 ) = (a 2 + 6 2 )*, so that 

L 2 = \og(a 2 + b 2 ). (12) 




where 

s T = log(l + b T /a). 

The s T are i.i.d., so the sum is normal for large t with 
mean t/j, s and variance to 1 , where \i s and a 2 are the mean 
and variance of the s T , by the central limit theorem. The 
system is thus log-normally distributed in the asymptotic 
limit and its moments are given by 



(x p t ) = x p a pt e pt ^ +pH<T ^ 2 . 



(13) 



Since we know that the first moment (xt) = x^a 1 is inde- 
pendent of the noise, we can conclude that (i s — — cr 2 /2 
and we have 



hi ~ We """""" ( 14 ) 
in the large t limit. Thus the Lyapunov exponents are 
given by 

L p =ploga- Hsp{p- 1), 

or 

L p = L° p - fi sP (p - 1) (15) 

where L p is the Lyapunov exponent for the unperturbed 
system. Notice that fx s < because the log function weights 
the negative values of bt /a more heavily than the positive 
ones. 

Expanding the log in the expression for /i s = (log(l + 
b T /a)) we find: 



(b T /a) 



k 



The (b T /a) k term in the expansion must be 0(b/a) k or 
smaller since b T /a can never exceed 1. Thus 

-p(p-l)(b/a) 2 /2- 



Lp — L p 



o(b/ a y 



(16) 



The error is 0(&/a) 4 if the noise is symmetric. In partic- 
ular, for the second moment, 



Ll + (b/a) 2 + 0(b/a) 3 



(17) 



in agreement with the exact value to second order. 

The approximation and exact results for a scalar sys- 
tem are compared to simulation in figure El below. This 
system converges in mean but has diverging moments for 
p > 3. The parameters for the simulation are a = 0.97 
and normal noise with b 2 = 0.05. The solid lines show the 
average of the moments (x P ) for p — 1,2,3 and 4, over 
10 6 runs. The dashed lines are the exact prediction pi(l 
and are shown only for p = 3 and 4. The crosses are the 
approximation l|l(j|) ; the inaccuracy for p = 4 is due to 
the expansion of the log. The initial value was xq = 1 and 
noises larger than a were not allowed. 



2.2.2 Approximate asymptotic distribution 

In this subsection we assume small noise, that is, |6 T /a| < 
1. This allows us to take logs and ensures that the mo- 
ments of b T /a are well behaved. We have 

log — = t log a + s T 



2.3 A = limit of multivariate case 

Returning to the multivariate system, we consider as a 
first treatment the A = limit. The system's evolution in 
this limit is given by 

x< = [nrr]x°- (is) 

T=l 
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time 

Fig. 2. Moment evolution in scalar system 



If the noises did not all have the same variance, the re- 
sult would be identical with b 2 replaced the average vari- 
ance 



2.3.2 Correlated noises 



When all the noises are correlated with the same variance 
(T noise) the calculation is similar except that both the 
sum on the {j T } and the sum on the {k T } in equation (|20|l 
give a factor of n. We thus obtain 



(|xf) = (nV)Vl 1 



(22) 



The expected value of x is 0. The expression for moments 
contains two or more occurrences of each B T ; the difficulty 
in its evaluation, and in general the difficulty of any multi- 
variate system, is that the noise matrices do not commute. 
However, when A = and the noise is white, the sum may 
be evaluated explicitly. Its value depends on the type of 
noise considered. In this section we show the details of 
how to evaluate such a sum; in later sections such steps 
will be skipped. 

In principle any moment of x could be calculated ex- 
actly, given a particular distribution for the noise ele- 
ments. Here we restrict our calculation to the second mo- 
ment for clarity and simplicity. 



2.3.1 Independent noises 

Wc first consider the case where all the noise elements vary 
independently The matrices B T do not commute so we 
must consider the full term by term expansion to evaluate 
the second moment: 



i xt i 2 )-E E E (Bk*B. 

1 ji32---jt kik 2 —kt 



2 

hh 



(19) 



r* Jr 1 r 2 r* ^.o\ 



where the expected value goes inside the sum because it 
is a linear function. All the elements of every B are inde- 
pendent, and we get a <Jj T fc T for every r = 1 . . .t when we 
sum on the fc's. This gives 



We note that this result can be also found immediately 
by transforming to the eigenspace, since in this case B is 
proportional to the matrix G of all ones which has A = n 
and all other eigenvalues equal to 0. 

The case of noises where only certain elements are cor- 
related provides an intermediate case between indepen- 
dent and correlated noises. The expressions are complex 
and are left for future work. 




Fig. 3. Second moment evolution for A = 



Second moment evolution in the noise-only case is shown 



i x ti2\ = u B ] ) 2 )!(B 2 ) 2 ) • . . . • ((B* ) 2 (x ) 2 ) in figure The simulations show the value of |x| 2 aver 



(20) 

When all the noise elements have the same variance b 2 , 
each of the t sums on j'2 , ■ ■ ■ , jt and i simply gives a factor 
of nb 2 . The remaining sum is just the norm squared of x°, 
and we obtain 



aged over 1,000 runs, in the independent noise case, and 
100,000 runs for correlated noise 2 . Here n = 3 and the el- 
ements of B were chosen from a normal distribution with 
variance 0.25. The simulations are compared to the pre- 
dictions of ||2U and (22| • 



,*|2\ 



(n& 2 ) 4 |x°| 2 



(21) 



2 A large number of runs is necessary in the correlated noise 
case because of the divergent moments 4 
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3 Properties of multivariate stochastic 
systems 



3.1 Log-normal character of distribution 



As we saw in H2.2I scalar stochastic systems with station- 
ary multiplicative noise are log-normally distributed with 
parameters proportional to time, so the system moments 
evolve as exp{tp(fi + pa 2 /2)]. While \i is typically nega- 
tive, for large p the positive pa 2 /2 term dominates and 
causes divergence. The effect of the multiplicative noise is 
thus to cause the system's pth moments to diverge for all 
p greater than some pq. 

While the components of multidimensional stochastic 
systems with multiplicative noise do not have an exact log- 
normal distribution, they retain the general log-normal 
character including the heavy tail and divergent moments. 
To be exact, any element of a product of t stationary ran- 
dom matrices is asymptotically log-normally distributed 
with parameters proportional to t |33II34| . Components of 
a multivariate stochastic difference system are thus linear 
combinations of log- normal variables with parameters pro- 
portional to t. Just as in the scalar case, therefore, multi- 
plicative noise in multivariate systems causes the system's 
moments to diverge. 




-1 -0.5 0.5 1 

log xf - log < xf > 

Fig. 4. Approximately log-normal distribution of the system 
state 



In the particular case of small noise and simple domi- 
nant A, the distribution of the elements of a multivariate 
system is very close to log-normal. This is shown in the 
simulation of figure Q] which presents a histogram of the 
log of 100,000 instances of xf° in a positive system whose 
A is given in l|63(l . The data were normalized by the ex- 
pected value (xf°). The solid line is a Matlab normal fit 
with n = -.0389 and a = 0.2725. The dashed line is the 
prediction of H4.3l and has fi = —0.348 and a = 0.2639. 



3.2 Relevant properties of A 

3.2.1 Simple dominant eigenvalue 

In this paper we only consider systems with simple, domi- 
nant A. Geometrically, the effect of A repeatedly acting on 
a vector is to bring that vector into the direction of u and 
to multiply its length repeatedly by A. The behavior of un- 
perturbed multivariate systems with a simple, dominant 
A is thus equivalent to scalar systems in the asymptotic 
limit. 

The requirement that A be simple and dominant is 
met in all nonnegative systems of interest (appendix 0), 
so our treatment of nonnegative systems is comprehen- 
sive. Although many arbitrary systems meet this condi- 
tion as well, some do not and we do not attempt to treat 
these cases. We also neglect systems with defective (non- 
diagonalizable) A, which form a set of measure 0, because 
the nonzero elements of A are impossible to determine 
exactly in most applications. 

3.2.2 Condition of A 

The effect of noise on a multivariate system, from a ge- 
ometric perspective, is to perturb both the direction and 
length of the vector x. Noise as a small perturbation means 
that a given noise matrix does not swing the x far from 
the direction of u or multiply |x| by a factor far from A. In 
this regime, the dynamics are well approximated by the 
dynamics of a perturbed scalar system. 

The regime of small noise, for the multivariate systems, 
is determined not only by the size of the noise elements 
but also by the sensitivity of the system to perturbation. 
There exist matrices whose eigenvalues and eigenvectors 
are violently affected by even a small perturbation to the 
matrix elements 35 , 36 . For a perturbation treatment, we 
need to know how much the dominant eigenvalue A and 
its eigenvector u of the system are perturbed by a given 
level of noise. 

The response of A to noise is characterized by a quan- 
tity k(X) called the condition of A. When k(X) is large A is 
said to be ill-conditioned, meaning that its response to a 
system perturbation is large with respect to the perturba- 
tion. Even a small noise causes moment divergence in sys- 
tems with an ill-conditioned A. Conversely, when k(A) = 1, 
A is said to be perfectly conditioned; its response to a sys- 
tem perturbation is the smallest possible and is on the 
order of the size of the perturbation. In systems with a 
well-conditioned A, the perturbation approximation is ap- 
plicable to relatively large noises. 

The change in A due to a small noise matrix B (small 
in the sense that 
|B| = 5 <C 1) is given by 

SX m v • (Bu) (23) 

to first order in S. Taking norms, we obtain the expression 
for the condition of A in the case of normalized u: 

k(A) - |v|. (24) 
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It is clear that k > 1 by the Schwartz inequality. The 
sensitivity of u to noise may also be calculated to first 
order[3S] and depends on the condition of A. It also de- 
pends on the gaps X — Xi between the dominant eigenvalue 
and the others, and is therefore related to the accuracy of 
the approximation 



A p « A p uv J 



(25) 



obtained by neglecting all Af compared to A p . In the limit 
that A2 — > 0.1251 is exact for all p and the sensitivity of u 
is minimized. 

The level of noise which qualifies as a small perturba- 
tion must therefore depend on k(A) and A — | A2 1 , which it 
does, as we will show. The effect of a given level of noise is 
the smallest in well-behaved systems with a well-behaved 
A and small eigenvalue gap. 



Matrix type 


Eigenvalue gap 


Condition k(A) 




A-|A 2 | 


of A 


A = aG 





1 


small a\ 


small 


close to 1 


normal A 


? 


1 


large a\ 


possibly small 


possibly large 



Table 1. Properties of some types of A. Recall that G is the 
matrix of all l's (mean value approximation). 



It is difficult to generally characterize the condition of 
A and the eigenvalue gap in terms of more physical prop- 
erties of the matrix A. What we can say is summarized in 
Table H] 

Systems close to the mean value approximation (recall 
A « aG where G is a matrix of l's in the mean value 
approximation; H1.2|) are sure to be well-behaved; however, 
some systems far from the mean value approximation are 
also well-behaved, as shown in figure El below as well as 
figures El an d E] 

The correlation between eigenvalue gap and condition 
number of A are demonstrated in the scatter plots of fig- 
ureOwhich show the eigenvalue gap versus k(A) for 10000 
randomly generated 5x5 matrices. The matrices were gen- 
erated from a normal distribution (top left), uniform dis- 
tribution (top right), uniform distribution with probabil- 
ity 1/2 and with probability 1/2 (bottom left) and uni- 
form distribution with mean 0.2 and variance 0.02 (bot- 
tom right). For the arbitrary matrices, only those with 
real A were accepted and the entries were normalized so 
that A = 1. Note the difference in the regions plotted. As 
shown in the figure, the likelihood that a given system is 
well-behaved is larger for nonnegative matrices than for 
arbitrary matrices, and larger still for positive matrices. 
See appendix ICl for further discussion. 



3.2.3 Limits on average element size 

Finally we note that in the case of nonnegative matrices, it 
is impossible to have a small A if the elements of A are too 



large. Many quite accurate bounds on the largest eigen- 
value of nonnegative matrices exist (see [HZj for a list); a 
relatively inaccurate but analytically tractable bound is 
the row sum bound, min^^^ Ay) < A < maXjQ^j Ay). 
This estimate implies that on average we need to take 



a < l/n 



(26) 



to keep A < 1 and ensure that the system converges in 
mean. This is exactly the asymptotic n result of |38] . and 
the result we would obtain in the mean value approxima- 
tion A « aG. 



3.3 Types of noise for multivariate systems 

For multivariate systems many different forms of noise are 
possible, distinguished by whether the elements are corre- 
lated and how large their relative variances are. In this pa- 
per we consider five cases which are analytically tractable 
and have some relevance to physical systems. The corre- 
lation rules for these cases are shown in table which 
provides a summary. 

For the correlation we consider three cases. Uncorre- 
lated noise means that the elements of the noise matrix 
vary independently. Totally correlated noise means that 
all the noise elements vary in the same way at each time 
step. For symmetric systems, we consider symmetrically 
correlated noise. 

For the variance we consider two possibilities. For ho- 
mogeneous noise, the variance of every element is identical 
and equal to b 2 . For proportional noise, the standard devi- 
ation of Bij is proportional to Ay by some factor q which 
we will take to be less than 1. 



4 Small noise as a perturbation 

In this section we determine approximate expressions for 
the moment Lyapunov exponents for multivariate systems 
subject to small noise using a perturbation treatment. We 
examine the dependence of the Lyapunov exponents on 
system properties, and discuss the accuracy of the ap- 
proximation. 

First let us reexpress the matrix product in JSJ: 



JJ(A + B T ) 

r=l 

E 



(27) 



E E 



Y*y* 



•Y 1 



Y'=A,B f Y t_1 =A,B t_1 Y 1 =A,B 1 

x 1 



E 

each y t =a.b t 



Y*y* _1 • • • Y 1 "' 



(28) 



meaning that each Y r in the sum can be either A or B r , 
for t — 1 . . . t. There are 2* terms in the sum; each term 
is a vector. 
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l^l-l^l 



Primitive 
matrices 




c< 1.4 

2 







Positive 






matrices 





















3 1.02 



Mean value 
approximation 




^ - l^l \ ~ \^\ 

Fig. 5. Eigenvalue gap versus condition of A in randomly generated matrices. 



Noise type 


Correlation rule 


Uncorrelated homogeneous (UH) 
Symmetrically correlated homogeneous (SH) 
Totally correlated (T) 
Uncorrelated proportional (UP) 
Symmetrically correlated proportional (SP) 


(BijBi'ji) =bi i 8 ii ,6jj, 
{BijBiij,) = b 2 (Sii>5jj> +S ij ,5 i ,j) 

{BijBi,,,) = q 2 (A ij ) 2 S ii ,S jj , 
(BijBiijt) = q 2 ( y A iJ ) 2 (8 ii ,8 JJ , +8iji8 Vj ) 



Table 2. Correlation rules for types of noise considered in this paper 



4.1 Perturbation expansion 

The perturbation expansion consists of considering only 
terms in l|28|l which have very few B's. For small noise, 
these terms make the only important contribution to the 
sum. Let us assume that this is so without justification, 
even before we define small noise. 

The reason that this strategy simplifies the calculation 
is as follows. Consider the evolution in time of the length 
and direction of a single term of l|28|l with few B's. In 
the asymptotic limit, a typical term with few B's has long 
strings of consecutive A's broken by single occurrences 
of B's. As far as the direction of such a term, the long 
strings of A's act to bring it parallel to u as previously 
mentioned (see (J2EJ0- When a B T acts on the term, the 
term lies almost parallel to u; even though the noise causes 
the term to point away from u, the next string of A's 



brings it back to the direction of u before another noise 
term occurs. The action of B r is thus independent of r. 
As to the length, a string of p A's simply multiplies the 
term length by A p ; and the B's multiply the length by 
some stationary random variable. 

In a term with few B's, therefore, the position of the 
matrices in the sum (|28|) is unrelated to their net effect on 
the term. Thus the matrices in the sum can be replaced by 
scalars, and the matrix product l|27|l becomes a product 
of scalars. To illustrate this, consider a typical term for 
t = 10 with a B r only in the r = 6 spot: 



AAAAB 6 AAAAAx° w (A 4 uv T )(B 6 )(A 5 uv T )x° 
w A 10 e 6 u(v ■ x°) 
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where we define the random variable 



In general, a term of the sum l|28(l that has long strings 
of A's and rn isolated B's {B Tl , . . . , B Tm } points in the 
direction of u and has length X t (e Tl ■ . . . ■ e Tm )u(v • x°). 
Such terms dominate the sum (Tggj) (see ET5j) and so the 
system state is given approximately by 

t 

x^u^-x )^ J](l + e T ). (30) 

T = l 

The random variables e T are i.i.d. and satisfy (e r ) = 0; the 
moments depend on the form of the noise. Notice that the 
numerator of e T is exactly equal to the first order change in 
A due to a small perturbation to A if^l and thus closely 
related to the condition n(X) J53J. The eigenvalue gap 
and thus the sensitivity of u is implicitly involved in this 
expression from the application of (|25|l . 

4.2 Criterion for small noise 

The simplest small noise criterion is 

(e P T ) « 1 

for all p. This is a rather complicated condition since the 
calculation of all the moments can be difficult for some 
forms of noise. Instead we choose a more restrictive (triple) 
condition, 

P(|e T | > 1) = 

(eP) ~ e*', p'<p (31) 
e 2 < 1 (32) 

Note that this requirement is not trivial as in the scalar 
case because the condition of A can be large. Less restric- 
tive conditions are possible but this will enable us to better 
understand the dynamics by taking logs and expanding in 
a power series in e. 

4.3 Moment evolution 

Using the perturbation expansion of section we may 
now present approximate expressions for the moments of 
a multivariate stochastic system. We do so by calculating 
the approximate Lyapunov exponents, proceeding from 
l|30(l exactly as in the scalar case of W2. 2. 21 with e T playing 
the role of b T and A the multidimensional analog of a. 
We thus find 

L p raplogA-Kp- l)<ln(l+e T )) 

^L° p +p(p-l) £ ^+0(e 3 ) (33) 

where 

L°=plog\ 



is the Lyapunov exponent for the unperturbed system and 
the error is 0(e 4 ) if the noise is symmetric. The system 
moments are 

(|X*| P ) M | V -X°|V L r 

« \(^)\^ 2PJS ^ 1+0 ^. (34) 

In particular, 

L 2 « L° 2 + e 2 (35) 
to second degree in e, and 

(|x*| 2 )«|<x*)| 2 e te2 . (36) 

Notice that to this level of approximation, first moment 
(norm) convergence is not distinguishable from conver- 
gence in mean. 



Noise Type 


e 2 


UH 

SH 
T 
UP 
SP 


v b 
2b 2 

*T 

^-Z)ij-(^tf) 3 K ?u ? +v i v j u i u j ] 



Table 3. Values of e 2 for the types of noise defined in 



To proceed beyond these expressions we must evalu- 
ate e 2 = (( v jf" ) 2 ), which we cannot do without specify- 
ing the form of the noise. The values of e 2 for the noises 
described in t|3.3l are easily calculated and presented in 
table |3 Here UP is uncorrelated proportional noise, UH 
is uncorrelated homogeneous, SP is symmetrically corre- 
lated proportional, SH is symmetrically correlated homo- 
geneous, and T is and totally correlated. In the propor- 
tional noise, / is a factor which depends on the distri- 
bution chosen; for example, / = 1 for normal noise and 
/ = 1/3 for uniform noise. In the case of symmetrically 
correlated noise, a symmetric A is assumed. 

The accuracy of the above approximations for the mo- 
ment evolution is demonstrated in figure |SJ In this figure, 
the log of the 2nd and 3rd moments of two randomly gen- 
erated systems are shown. The plots are normalized by 
the expected value (unperturbed value) of the system and 
show only the noise part. The solid line is the average over 
10,000 runs of the simulation, and the dotted line shows 
the analytic prediction of l|34|) . At left is the positive sys- 
tem of (|B3|I subjected to uniform UP noise with q = 0.5. 
Note that the asymptotic limit is reached almost immedi- 
ately in this system. At right is an arbitrary system with 
simple dominant A subjected to normal UH noise with 
b = 0.1. This system has large transient behavior before it 
settles in to its asymptotic limit around t = 20. The an- 
alytic prediction, which cannot account for the transient, 
has been artificially placed to demonstrate the asymptotic 
accuracy of the slope. 
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Noise Type 


L 2 


UH 


L 2 f 






SH: 


L 2 f 




+ A— 


T: 


L 2 f 




+ q 2 



Table 4. Approximate value of L 2 for three types of homoge- 



neous noise 



4.4 Dependence on system size 

We can now explore the n dependence of the moments. Be- 
cause the small noise case is important for applications, 
we present a detailed discussion of the size dependence 
based only on the expressions developed thus far. A dif- 
ferent discussion of the n dependence for larger noises is 
presented in WiA\ For simplicity, we consider only the sec- 
ond moment in this section. 

As we show, independently varying noises "interfere" 
with each other and diminish the effect of the noise, com- 
pared to the unperturbed system. Thus, the effect of the 
noise decreases as n increases in the case of uncorrelated 
noise. There is no n dependence to second order, how- 
ever, in the case of totally correlated noise. Symmetrically 
correlated noise provides an intermediate case. 

We also show that for noise proportional to the system 
elements, as the system deviates from the mean value ap- 
proximation and in particular becomes closer to diagonal, 
the destructive interference is decreased and the noise has 
a greater effect. 



4.4.1 Mean value approximation 

As a first simplification, we consider the mean value ap- 
proximation where A « aG. In this case, Vi ~ Ui m l/\/n 
for all i, and A w na. We consider homogeneous noise (all 
the noise elements have the same variance, which is al- 
most equivalent to proportional noise in the mean value 
approximation) with variance b 2 = q 2 a 2 , q < 1. 

Using (|5o")l and table the values of L 2 for three types 
of homogeneous noise are easily computed in the mean 



value approximation and are shown in tabled These ex- 
pressions are to be compared to the scalar case L2 ~ 
L\ + q 2 (equation [T7J) . 

Note in particular how the noise effect (the q 2 term) 
is divided by a factor related to the number of indepen- 
dent elements of the noise. This destructive interference is 
not surprising when we consider why multiplicative noise 
processes generate the anomalously large events which 
make up the heavy tail of the log-normal distribution. 
The anomalous events result from a long sequence of large, 
positive noises |39j . When there are n 2 independent noises 
per time step, as opposed to 1, anomalous events are rarer. 
However, when all the elements of noise vary identically, 
the effect of the noise is the same in scalar and multidi- 
mensional systems. Note that symmetric noise provides 
an intermediate calculable case; there are n(n + l)/2 in- 
dependent components in a symmetric noise. 

The dependence of the noise effect on the number of 
independent elements of the noise is demonstrated in the 
simulation of figure [7| By tableland equation 



(N 2 )/|(x)| : 



Sfq 2 



in the asymptotic limit, where / = l/n 2 for UH noise (all 
noise elements independent) and / = 1 for T noise (all 
noise elements correlated). Accordingly we plot 
<7 -2 (log|x| 2 — log|(x)| 2 ), averaged over 100,000 runs, for 
various values of n. On the left is the plot for UH noise, 
where the lines should have slope l/n 2 ; on the right is 
the plot for T noise where the slope should be 1 for any 
n. The agreement is excellent. The A's for these systems 
were randomly generated from uniform distributions with 
small variance. The initial state was a vector of l's, and 
q = 1/4. Because we are in the mean value approximation, 
A2 is very small for these systems and the asymptotic limit 
for t begins almost immediately. 



4.4.2 Deviation from the mean value approximation 

We examine the effect of a deviation from the mean value 
approximation on Li in the case of uncorrelated propor- 
tional noise. The result is that the noise effect roughly 
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increases the larger the deviation, as the 1/n 2 damping 
caused by the independent noise elements is mitigated. 
This is because we have assumed that the typical size of a 
noise element is proportional to the corresponding element 
of A (UP noise), so that small entries of A contribute little 
to the interference effect of the independent noises. 

In this subsection we assume that the approximation 
is accurate for p = 2 so that A?, w X 2 u 2 v 2 . Recall 
that (|25|l is generally more accurate the closer A is to the 
mean value approximation, but it can be accurate even if 
the variance of the Aij is large, as discussed above. 

With the above approximation we have 

L-2 « In A 2 



2 4 

q w 



(37) 



system does not immediately reach its asymptotic limit; 



it has | Aa | 
system. 



A/2 as opposed to |A 2 | < 1 for the MVA 



— mean value system, w 
other system, w 2 = 0.3 



1/n = 0.2 



where we define 



Comparing (|37() to the mean value case, we see that the 
1/n factor is replaced by w 2 . This quantity satisfies 

1/n < w 2 < 1 



time 

Fig. 8. Larger noise effect of a deviation from the mean value 
approximation 



since v • u = 1. The lower bound is achieved in the mean 
value case; the upper bound is achieved when A is diago- 
nal. w 2 is thus a rough measure of the deviation of the 
from the mean value approximation; it generally increases 
as the variance of the increases. 

Larger values of w 2 for deviations from the mean value 
approximation thus mean a larger noise effect. This ef- 
fect is demonstrated in figure |H1 which is analogous to 
figure and compares the noise parts of the second mo- 
ment q ~ 2 (log | x | 2 — log | (x) | 2 ) of a mean value approxima- 
tion (MVA) system and of a system with the same A and 
n = 5 but whose elements have a much larger standard 
deviation. UP noise was considered, so the slope should 
be w 4 by {33J. The dashed line system has w 2 = 0.2888, 
w 4 = 0.0834 while the MVA system has w 2 = 0.2002 « 
1/n and w 4 = 0.04; a rough linear interpolation fit to the 
data shows an asymptotic slope 0.095 for the dashed sys- 
tem and 0.0385 for the MVA system. Note how the dashed 



4.4.3 Large n limit and homogeneous noise 

When the noise is homogeneous we can apply results on 
spectral theory of matrices to study the n dependence in 
the large n limit without appealing to a mean value ap- 
proximation. In the case of a symmetric matrix with en- 
tries drawn from a distribution with mean a and variance 

X = na + a 2 /a 
on average |40j . For an arbitrary (asymmetric) matrix |38| 

A ~ na, 

which is really just the mean value approximation. 

We thus obtain table for the n dependence. Recall 
the notation b 2 for the variance of the homogeneous noise. 
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Arbitrary system L2 ~ L% + b ^ 
~ 11 2 — 2 

Symmetric system L 2 ~ Dj + (n . +2 " n ^ /a . )/2 

Table 5. Large n dependence of arbitrary and symmetric sys- 
tems with homogeneous noise. 



Note again the n 2 damping in the arbitrary system, and a 
damping on the order of n(n + 1)/2 in the symmetric case, 
in agreement with the previous analysis. We have assumed 
independently varying noise in the arbitrary system, and 
symmetrically varying noise in the symmetric system. 

4.5 Approximation justification, accuracy, failure 

The justification for equation l|30|l in the small noise ap- 
proximation is as follows. Expand the product T (1 + s T ) 
into a sum. The typical size of the random variable e T is s, 
and the largest contribution to the sum comes from terms 
with /c max ev's, where 

is the binomial expected value. The brackets denote the 
closest integer. This means that the largest terms in the 
sum come from terms of (|28[1 with fc max B's. From symme- 
try considerations it is clear that in the asymptotic limit, 
the average separation d between two B's in a term with 

fcmax B S is 

Id) = 

w 1 + te' 

which is large for small e and asymptotically independent 
of t. Furthermore, in a term with fc max w te B's, the sep- 
aration satisfies 3 

P{d < d ) = 1 - (1 - - ed 

in the asymptotic limit, which is small for small e and 
independent of t. Therefore, the important terms of l|28[l 
for small e are those with a few B's separated by long 
strings of A's for all t 3 . 

However, this analysis does not tell the entire story. 
The accuracy of the perturbation approximation is in fact 
much higher than one would expect from the above cal- 
culation. To understand this, consider a term of the sum 
l|28l) with many B's. This term's direction is impossible 
to determine in general because each noise matrix trans- 
forms it arbitrarily. There are many such terms and they 
are all affected by a different set of noise matrices. Their 
directions are thus widely distributed in R™ and mostly 
cancel out in the sum. 

When e is not small terms in the sum (|28[) with many 
B's become important. This causes the perturbation ap- 
proximation to be inaccurate for two different reasons. 

3 Note that simulation of divergent moments in a convergent 
system for large t may not seem accurate, because as t increases 
the probability of an anomalous event becomes very small. A 
very large sample space is necessary to obtain an accurate sim- 
ulation for large t; see the discussion in 4 



First, when B's are adjacent, the approximation of replac- 
ing B by e T is poor; second, when there are many strings 
of only a few adjacent A's, both replacing A by A 2 and 
B r by e T can be inaccurate. The relative importance of 
these two inaccuracies can be different. For example, the 
accuracy of the A factor is independent of n while the 
accuracy of the B factor decreases as n increases. 

It is difficult to determine a cut-off where e becomes 
large. The overall error may be much smaller than the 
error of each term of the sum (|28(l , because the deviations 
of the terms may lie in different directions and cancel out 
in the sum. It is clear that the cut-off depends on how 
quickly A p brings a random vector into alignment with u, 
but even this is a complicated function of the eigenvalue 
gap and the condition of A ■ To account for large e and 
handle the contribution from neighboring B's accurately 
for large n, we develop a different approximation in the 
next section. 



5 Arbitrary noise using iteration 
approximation 

We now present a different method, the iteration tech- 
nique, which can be applied to find an approximate value 
for the second moment for small or large noise in well- or 
ill-conditioned systems. 

For homogeneous noise (that is, all the noise elements 
having the same variance), the approximation can be ex- 
tended to any level of accuracy for any noise. Unfortu- 
nately, for other forms of noise including proportional noise, 
only the first approximation is applicable. 

In addition to providing a way to treat systems where 
the noise effect is not small, this technique is able to detect 
the explicit n-dependence of the noise effect. This effect 
is very slight in the small noise case, but quite important 
for larger noises. 

The general strategy of the method is to express (|x t+1 | 2 ) 
as a time-independent function of {(|x*| 2 ), (|x* _1 | 2 ), . . .} 
in the asymptotic limit. A similar technique was indepen- 
dently developed in |41| for other applications. 



5.1 First approximation 

The first approximation of this method consists of apply- 
ing the relation 

^(^■) 2 «A 2 V, (38) 

where u has been normalized to have length 1, for all r, 
even r = 1. In this approximation we can express (|x* +1 1 2 ) 
as a ^-independent function of | 2 ) alone, as we will see. 
We thus define x^ = Ax' -1 and = Bx* _1 , so that 
x* = x^ + x^ . We have 

(|x*| 2 ) = (|x^| 2 ) + (|x* 3 | 2 ); (39) 
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the cross term is zero in expectation because there is one 
power of B* . We will establish a matrix recurrence relation 



and the second moment evolves as fj,\, where the largest 
eigenvalue \i\ of is 



A I / 
J+l|2\ 



(40) 



where the elements of Mi (subscript 1 for first approxima- 
tion) are independent of time. The asymptotic behavior of 
the second moment is (|x*| 2 ) ~ /i* , where \i\ is the largest 
eigenvalue of Mi, and 

L 2 « ln/ii 

for the Lyapunov exponent. 

Using the new notation on the recurrence relation, we 
have: 

<|x t+ T) = £ ((A + B*+% {x\ + x^ (41) 
in' 

(A + B t+1 ) if (x A +x B ) f ) 



or 



|x t+1 | 2 ) = (lAx^l 2 ) + J^A ij A ij/ ((x B ) j (x%) r ) + (42) 
in' 



because the noise is white with mean 0. This is the sim- 
plest form we can obtain without considering particular 
types of noise. 



5.1.1 Homogeneous noise 

Recall that homogeneous noise is a type of noise in which 
all the elements of B have the same variance. In this case, 
equation 1|42[) becomes 

<|x t+1 | 2 ) « A 2 (|x^| 2 ) + (X 2 f v /n){\ X %\ 2 ) + (43) 
+ (n k b 2 f u )(\^ A \ 2 ) + (n k b 2 )(\^ B \ 2 )) 

where we introduce the notation 

v 2 , UH noise 



fv 



fu 



as well as the factor 



(EiVi) 2 , T noise 

1, UH noise 
(UiU t ) 2 , T noise 



1, UH noise 

2, T noise 



(44) 
(45) 

(46) 



to account for the difference between independent (UH) 
noise and correlated (T) noise. Multiple steps have been 
skipped in obtaining equation l|43|l . including the use of 
(|25() with p = 2 on the first term and l|38|) on the others. 
We thus obtain 



M" H = 



A 2 f v X 2 /n 



nf u b n k b 



kl,2 



X 2 + n k b 2 + J(X 2 -n k b 2 ) 2 + 4X 2 v 2 b 2 , . 

Ml = ^ " ■ (47) 

Recalling that for UH noise v 2 b 2 /X 2 = e 2 , while for T 
noise b 2 (E i u i ) 2 (E i v l ) 2 /A 2 = e 2 , we thus have 



T UH , rx2 / 1 + n k b 2 /X 2 + - n k b 2 /X 2 ) 2 + Ae 2 



ln[A 2 (- 

L° 2 + e 2 + eH^-h + 0(s% 



)] 
(48) 



where L\ = 2 log A, and the approximation in the second 
line is valid in the limit of small e 2 and small nb 2 / X 2 . The 
main difference between this expression and the pertur- 
bation expansion is that we have taken into account the 
effect of two neighboring B's, which produces a factor of 
n. The n dependence enters only in the second and higher 
order terms; this expression agrees with the perturbation 
approximation 135fl to first order. 



5.1.2 Proportional noise 

In the case where the noise elements satisfy bij = qAij 
with q < 1, we apply 1351 with p = 1 and proceed as above 
to find 



L 



UP 



ln[A 2 ( 



1 + q 2 w 2 + - q 2 w 2 ) 2 + Aq 2 w 4 



)] 



1 



2\6 



L\ + {qw 2 ) 2 + {qw 2 )\— - -) + 0(qw 2 ) 



w 

where w 2 was defined previously f M4.4fl as ^2jV 2 u 2 , and 
the approximation in the second line is valid in the limit of 
small qw 2 . As expected this agrees with the perturbation 
approximation result (|35J) for proportional noise to first 
order. 



5.2 Further approximation 

The above treatment is completely accurate in the way 
it handles the B for homogeneous noise. Any inaccuracy 
stems from using the approximation A r ps A r uv T on the 
A for r = 1. We can improve on this inaccuracy to any 
desired degree, as explained below. Unfortunately, any ap- 
proximation past the first order is only applicable to ho- 
mogeneous noise (all variances the same) and not to pro- 
portional noise or any other form with different variances. 
For the remainder of this section, therefore, only homoge- 
neous noise will be considered. 



5.2.1 Second approximation 

To illustrate the idea, we begin with a second approxi- 
mation wherein l|38|) is assumed to be accurate for r = 2 



14 



Dennis M. Wilkinson: Moment instabilities in multidimensional systems with noise 



and higher, but not r = 1. In this second approximation, 
A's which occur "alone" (surrounded by two B's) in an 
element contribute a factor aiA 2 instead of just A 2 . 



We now break x into x 



^AB 



analo- 



gously to l(39|) . where x^ are the terms beginning with 
AA, etc. Proceeding just as above, we find that 



with 



(Ix^l 2 ) 

(ix^ +i i 2 >; 


«M 2 


(Ka 

(\*AB 

V <b4l 


( A2 


A 2 /«r 





= 





f v ai\ 2 


\nfub 2 


nf u b 2 


n k b 2 



where f u , f v and k were previously defined in 1(45(1 . 1(44(1 
and 1(46(1 and account for the difference between UH noise 
and T noise. 

The second moment will diverge when the largest eigen- 
value of Ma is greater than 1. This eigenvalue is the largest 
root of the equation 

V? - M 2 (A 2 + n k b 2 ) + ^\\b 2 {n k - f u f v a{) + 
b 2 \if u f v (l~ ai )=0. 

Notice that in the limit ol\ = 1, that is, the limit that the 
first approximation is accurate, we recover the character- 
istic equation for the first approximation 147(1 . 



5.2.2 Higher order approximation 

We can extend the above procedure to any level of accu- 
racy. Define a vector a by 



1 



fufv^i r 



ab 



The elements of a are the successive corrections to (138(1 . 
As r increases, a r tends to 1 because A^ 3> Af becomes 
very accurate for large p. The Lyapunov exponent of the 
system is given by the log of the largest eigenvalue of 



M r 



( A? Afi 
A? 

- 1 rv„ 





\nf u b 2 nf u b 2 nf u b 2 



\ 2 22 



nf u b 2 







Ai«i/„/n 
n k b 2 



(49) 



in the large p limit. The characteristic equation for this 
matrix can be expressed iteratively as in ^15.31 below, but 
the largest eigenvalue must be computed numerically. This 
method is exact for any noise and any A with a simple, 
dominant eigenvalue, however ill-conditioned A may be. 

The accuracy of the higher-order approximation method 
is demonstrated in the simulation of figure In this fig- 
ure, the second moment of a system with a very poorly- 
behaved A feauation l64(l . subject to normal UH noise with 



simulation 
prediction 




time 

Fig. 9. Accuracy of approximation for very ill-conditioned sys- 
tem with large noise 



various b 2 , is simulated. The solid lines are the average 
over 10,000 runs of the simulation, and the dotted lines 
are the analytical prediction of H5. 2. 21 for r — 6. The A in 
this system has v 2 — 170.51 and so virtually any noise is 
not treatable using the perturbation approximation. The 
slopes of the dotted lines are 0.255, 0.833, and 1.794 for 
b 2 = 0.04,0.25, and 1 respectively. Compare to the per- 
turbation approximation which estimates 2.02, 3.85 and 
5.24. This system is convergent in mean with A = 0.95; 
notice how quickly the moments diverge even for small 
noise because v 2 is large. 



5.3 Large noise limit 

Using the results of the higher order approximation, we 
can obtain an approximation for the Lyapunov exponent 
for second moment evolution in the large noise limit. 

Note that situations where a small noise has a large ef- 
fect because the system's dominant eigenvalue is ill-condi- 
tioned (0 are not treatable using the formalism of this 
section, for the reasons given below. 



5.3.1 Criterion for large noise 

A good estimate for the onset of the large noise regime 
can be obtained by comparing the largest eigenvalue of A 
to the average largest eigenvalue of the B. 

For independent (UH) noises with variance b 2 chosen 
from a normal distribution, the magnitude of the largest 
eigenvalue of B is given on average by |4"2~] 



Ab — b\/n, independent noises 



(50) 



For correlated (T) noises chosen from a normal distribu- 
tion, the matrix B is simply a normal random multiple of 
the matrix G of all ones. G has largest eigenvalue n and 
so the largest eigenvalue of B is on average 



Ab = bn, correlated noises 



(51) 
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The large noise case corresponds to As ^> A, that is, 

n k b 2 > A 2 (52) 

where k = 1 for independent (UH) noise and k = 2 for 
correlated (T) noises. 

5.3.2 Lyapunov exponent 

To find the Lyapunov exponent for second moment evolu- 
tion in the large noise limit that n k b 2 ^> A 2 , we introduce 
the small parameter 



5 = Sis ^ L 
no 



(53) 



We will appeal to the iteration treatment of M5. 2. 21 which 
was accurate for any noise. Recall that the Lyapunov ex- 
ponent is given in this treatment by the log of the largest 
eigenvalue fi r of the matrix M r ( equation I49f) . where r 
is any integer. The approximation is more accurate the 
larger r is, but as we will see, in the large noise limit there 
is no need to consider large r. 

Note that the parameters {ai}, f u , and /„ enter into 
the calculation of \i r (equation I49fl . If these parameters 
are large, they can ruin our expansion since they multiply 
S. We therefore assume that they are O(l). Note that this 
amounts to assuming the system is well-behaved; this is 
why, as noted above, this expansion is not applicable to 
ill-conditioned systems. 

The characteristic equation for M r+ i can be written 

= (A 2 - M )hxD r + a r X 2r b 2 f u f v ] + X 2r b 2 f u f v 

where fj, are the eigenvalues and the D r are defined recur- 
sively by 

D r +i = -^D r + a r X 2r b 2 f u f v 

with 

-fx f v ai\ 2 /n'' 
nf u b 2 n k b 2 — fi 

Keeping only terms to first order in S in the character- 
istic equation above, we find that the largest eigenvalue 
/! r +i is given approximately by 



Do = det 



fr;,2 



independent of r, so that 



L 2 *\ogn k b 2 + 6^^ + 0(6 2 ) 
^\ogn k b 2 + \ 2aifufv ' 



i 2k b 2 



0(A 4 ). 



(54) 

(55) 
(56) 



Note that the zeroth order term corresponds to that found 
in the A — limit by different means in ^2.31 

The range of applicability of the large noise approxima- 
tion is demonstrated in figure where the second moment 
is plotted against the prediction of H55fl . The A of i|tj3f) 
is used, which has n = 5 and A = 0.96. Thus b = 1 is 
well within the large noise regime, as shown in the figure. 
Even b — 0.25, which is not within the large noise regime, 
is reasonably well predicted by this approximation. 



simulation 






— large noise prediction 








b = 1 








b = .25 





5 10 15 20 25 30 

time 

Fig. 10. Accuracy (and inaccuracy) of large noise prediction 

6 Critical value and stability diagram 

As a general rule, large deviations from the average be- 
come reasonably likely when the noise is large enough that 
the second moment diverges. The onset of second moment 
divergence therefore marks a threshold between two types 
of behavior and defines a critical value of the size of the 
noise. 

In the case of a homogeneous noise, in which all the 
noise elements have the same variance b 2 , the critical value 
can be simply expressed as the value of the variance where 
the second moment Lyapunov exponent L 2 equals 1 . For a 
proportional noise, the critical value is the value q c of the 
constant of proportionality (see table EJ) for which L 2 = 1. 
The approaches used in this paper allow for a detailed 
treatment of the critical value b 2 in the case of homo- 
geneous noise, but unfortunately not for proportional or 
other forms of noise, which is left as a topic for future 
work. Proportional noise is only discussed in the mean 
value case where it takes the same form as homogeneous 
noise. 

Throughout this section we will assume that the sys- 
tem is well-behaved, that is, f u , f v and the {ai} are close 
to 1. Recall that f u and /„ account for the difference be- 
tween independent and correlated noises, and the {ccj} 
measure the accuracy of the approximation (|38|l for suc- 
cessive powers A 1 . 



6.1 A 2 -> 1 limit 

A limit of particular interest when considering the critical 
value is A 2 — > 1, where, as we will see, the critical value 
drops sharply to 0. Since only small noise is required to 
cause divergence in this limit, we can apply the first ap- 
proximation of the iteration treatment, in particular equa- 
tion 03 to find the critical value 



( & c)small 



1 



1-A 2 



(57) 



where k = 1 for independent noises and k — 2 for corre- 
lated noises. The sharp dropoff to of the critical value is 
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evident from this expression and demonstrated in fjgure lTD 
below. Note that when n = 1, in which case f u = f v = 1 
and A = a we retrieve the scalar result b\ = 1 — a 2 . 

Our above expression for 6 2 should coincide with that 
implied by the perturbation approximation of <21 Using 
equation [2U combined with table |3| and definitions i|45[) 
and Ij44(l for f u and f v , we obtain the perturbation ap- 
proximation result 



( & c)pert 



1-A 2 

fufv 



which is easily seen to be equivalent to l|57[) to first order 
in (1-A 2 ). 



6.2 A 



limit 



When A 2 is small and the system well-behaved, the second 
moment will only diverge if the noise is large. In this case 
the limit can be obtained from the large noise treatment 
of H5.3I in particular equation (|54")l whence we find 



2\ 

) large noise 



(58) 



For well behaved systems with a± ~ 1, this expression is 
almost equivalent to the small noise expression l|57(l ! The 
scalar result is again retrieved from this expression. We 
also note that when A = 0, we obtain the critical value 



(bl) 



A=0 



l/n\ 



which was obtained in a different way in section f2.3l 



6.3 Stability diagram 

In the case of a well-behaved system with k 1, the 
functional form of the critical value is the same for small 
and large values of A 2 . We thus propose the following ex- 
pression 



b.= 



1 



1-A 2 



(59) 



for the critical value for all ranges of homogeneous noise in 
well-behaved systems. Using this expression we can make 
a phase plot for the stability regions of the system, as 
shown in figure ITTI for n = 5 and independent (UH) noise. 

As for proportional noise, we can produce a stability 
diagram to compare with the homogeneous noise case in 
the mean value limit. This diagram indicates how large, 
as a fraction of the size of the unperturbed elements, the 
noise must be to cause divergence. 

Recall that the size of a proportional noise was defined 
by the factor q, the constant of proportionality between 
the typical noise size and the average element of A (see 
table HJ. In the mean value approximation, the largest 
eigenvalue is given approximately by A « an and we thus 
obtain the critical value 

q c = -rb c . 

This relation produces the phase plot of figure El f° r the 
constant of proportionality q c as a function of A. 




Fig. 11. Second moment stability diagram for n = 5 




Fig. 12. Stability diagram for proportional noise in mean value 
approximation, n — 5 



6.4 n dependence of critical value 

Expression (|59|l can be used to study the n dependence 
of the critical value. The n dependence is weak for large 
A, but it is strong when A is small and a large noise is 
required to cause divergence. 

In the limit A — > 1 , where only a small noise is needed 
to create second moment divergence, it is seen from ex- 
pression l(5^|l that the n dependence is quite weak. Indeed, 
the expansion (|48|l for small noise showed that n depen- 
dence enters only in the second order term in this limit. 
The effects of different forms of noise in this case are quite 
similar. 

The n dependence of the critical value for small, inde- 
pendent (UH) noises is demonstrated in figurcEl The fig- 
ure plots the value of the second moment at t = 30 for sys- 
tems with various 6, normalized and averaged over 300,000 
runs. Because of the normalization, the initial value (at 
t = 0) of the system's second moment was 1. The criti- 
cal value 6 2 is thus indicated for each n by x-coordinate 
for which the plotted curve's y-coordinate begins to ex- 
ceed 1. Compare the values of 6 2 given by the simulation 
to the analytic estimates b\ = 0.0365, 0.0326, 0.0288 for 
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Fig. 13. n dependence of critical value 



n = 3,6,10, respectively from equation (|5§|) . The values 
of the elements of the As used in this plot were generated 
randomly from a normal distribution 4 , and the matrices 
were normalized to have A = 0.98. 

In the opposite limit of A — > and large noise, the n 
dependence is quite strong, as shown in l|55(l . It is in this 
limit that the difference between independent and corre- 
lated noises is quite marked, due to the "destructive in- 
terference" phenomenon of independently varying noises 
discussed previously in H4.4I 



6.5 Comparison to convergence bounds 

The critical value (|59|l provides a much more accurate es- 
timate of the "safe" level of noise for which the second 
moment does not diverge than do the convergence bounds 
of appendix IdI We bring this point up because the tradi- 
tional mathematical approach to stochastic stability is to 
use bounds. 



These bounds are b 2 < 



for the second moment 



1721) and b 2 < ^ ^ for any moment (|71|) in the large 
n limit. For typical well-conditioned systems with v 2 rel- 
atively close to 1 (figure 03 appendix Q , it is clear that 
the critical value is much less restrictive than either of 
the bounds. That is, the bounds stipulate that we must 
take a very small noise to guarantee convergence of the 
second moment; but the critical value indicates that the 
second moment will converge for a much larger range of 
noise. When v 2 is large, A and thus the matrix A are ill- 
conditioned and the norm of A is typically much greater 
than A, so the above bounds are not accurate. 



7 Results are generally inapplicable to 
continuous limit 

As a last subject, we discuss the continuous limit of our 
stochastic system. The results of this paper are not gen- 

4 the matrix was not accepted if it did not have a simple 
dominant eigenvalue 



erally applicable in the continuous limit because there is 
no such thing as small noise, in the sense we have used, in 
the continuous limit. Of the cases we have considered, the 
discrete result is only applicable to the continuous limit 
in the mean value approximation that A ss aG. 

The reason that the relative size of the noise depends 
on the time scale is that the correct limit of a white noise 
process has standard deviation proportional to ^/di |43| . 
Thus, the noise necessarily dominates as dt — ► 0. To il- 
lustrate this point, consider a particle moving in a one 
dimensional diffusion process 

dx = ((a — \)dt + bdw)x 

where dw is a Wiener process. When we consider the sys- 
tem's average motion on a large time scale, the particle 
generally progresses along the curve xoe^ a_1 ^. However, 
on very small time scale, the motion is completely erratic 
because it is dominated by the noise. 

For multidimensional systems, the continuous limit of 
J2J is the stochastic differential equation (in the Ito sense) 



dx = [(A - l)dt + dB]x 



(60) 



where dB = bdW is a matrix of Wiener processes with 
mean and standard deviation proportional to b\fdi, and 
I is the identity matrix. Again, in the dt — > limit, the 
motion is completely dominated by the noise and the vec- 
tor x is transformed erratically around in R™ . The system 
can never become aligned with u because the large noise 
causes it to couple with the other modes of A. Only when 
A2 — ► does the system become aligned with u and behave 
similarly to the perturbation approximation, above. 

7.1 Correspondence between continuous and discrete 
results in the mean value approximation 

For correspondence between the discrete and continuous 
cases we consider a system in which A2 = 0: the mean 
value limit that A = aG. For this A and totally corre- 
lated noise, an analytic solution to H60[) is possible be- 
cause A and dW = Gdw commute [44U26] . Here dw is a 
one-dimensional Wiener process and A = na is the only 
nonzero eigenvalue of A. The solution to Ij60|l is 

X(t) = e (A-I-ft 2 G 2 /2) t +6G lUx(0) ^ 

where w = L dw is normal with variance t. From this it 
is straightforward to calculate that 



(|x(*)| 



-Gx(0) 

n 



in the asymptotic limit, and the moment Lyapunov expo- 
nent is 

TL 2 b 2 

l p = -p5 + p{p-l) — 

where we have taken A = na = 1 — S. This can be com- 
pared with the discrete result for the mean value limit and 
totally correlated noise: 

L p » - P S + p{p - 1)^(1 - 26) + 0{8 2 ) + 0(6 4 ), 
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where we have applied 



Ui 



" x / 2 . In the limit of 



small time step the expressions are equivalent to lowest or- 
der. This same analysis can also be performed for a scalar 
system where there are no other modes to couple to. 



7.2 Failure of discrete result in the continuous limit 

When there are nonzero modes for A other that A, the 
discrete result should not, and does not, correspond to 
the continuous limit. This can be verified by comparison 
to the result of |27| for small noise moment Lyapunov 
exponents of arbitrary two-dimensional linear stochastic 
differential equations. This result, for white noise, is 



l p = -pS+pj! — +p 2 72" 



+ pO(b 2 ) + 0(p 2 ) (61) 



where we take A = 1 — 5. The 7 factors depend on the 
form of noise considered. 72 depends only on the domi- 
nant eigenmode, while 71 depends on both eigenmodes. 
To proceed we assume UH noise for definiteness, wherein 



one can show that 72 
UH noise is 



The discrete result (|35f) for 



L,, 



- P S + p{p - 1) — (1 - 26) + 0(5 2 ) + 0{b A ). (62) 



as described in the text from distributions with low vari- 
ance, and it would be a waste of space to present them 
exactly (see discussion in H3.2fl . 

The matrix which was used to generate figure ^ an d 
many others as noted in the text is 



A = 



/ 0.1795 0.0861 0.1860 0.0924 0.1661 \ 
0.1429 0.1680 0.0517 0.2626 0.3272 
0.3558 0.0127 0.2797 0.0221 0.3227 
0.2766 0.2654 0.1611 0.0408 0.0745 

\ 0.3539 0.3059 0.0596 0.2933 0.3147 / 



(63) 



This matrix has largest eigenvalue A = 0.966, second largest 
eigenvalue | A2 1 = 0.228, and v 2 = 1.10 as computed by 
Matlab. This matrix is thus quite well-behaved as defined 
in ffHl 

The matrix with ill-conditioned A used to generate the 
plot of figure |5] is 



/ 0.5086 0.3496 
-0.6168 0.1553 
-0.5526 0.0069 
0.4805 0.8053 
\ -0.4307 0.8960 



0.0795 

0.5224 

0.0008 

-0.5502 

0.0255 



-0.2044 

-0.0293 

-0.3189 

0.6173 

0.1454 



-0.3530 \ 

0.0137 

0.4345 
-0.3041 

0.6965 / 



(64) 



with largest eigenvalue A = 0.950, second largest eigen- 



value 
lab. 



A 2 | =0. 



and v = 170.3 as computed by Mat- 



Comparing this expression to the continuous version (|61|) , 
we see that the 1 — 25 factor on the noise term accounts for 
the difference between discrete and continuous evolution, 
as in H7.ll and that the p 2 ^i\ +pO(b 2 ) in (f^Tfl probably 

corresponds to the p(p — 1)^- term in l)62[l. Note that 
the p(p— 1) form is present in both continuous scalar and 
T noise cases and is typical of log- normal distributions. 

However, the term in l)61|l proportional to 71 is com- 
pletely absent in the discrete result; moreover, it depends 
on A2 and its eigenvector which have no effect on the small 
noise discrete system. This term shows how the solution 
is coupled to all modes, not just the dominant one, in the 
continuous limit. In fact, for UH noise, one can show (see 
H67|)1 that 71 = 1 — v 2 ; in the mean value limit v 2 = 1 and 
the contribution of the second mode is 0. 
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A Matrices used to generate figures 

All of the unperturbed matrices A used to generate the 
figures of this paper were randomly generated. Two par- 
ticular As are presented here. The others were generated 



B Reduction of nonnegative stability analysis 
to primitive systems 

The reason that A is simple and dominant in all nonnega- 
tive systems of interest is that we need only consider sys- 
tems with primitive A, and primitive matrices have the 
above property by the Perron-Frobenius theorem. Stabil- 
ity analysis of any nonnegative system whose matrix is 
not primitive reduces to analysis of primitive subsystems. 

More precisely, nonnegative matrices which are not 
primitive may be either reducible or irreducible imprim- 
itive. Reducible matrices are those which can be written 
in the form 

where C and D are square, by renaming the indices 
Stability analysis reduces to analysis of the subsystems C 

and D, because ( ^ ^ ~\ — ( ^ ~\ . A similar reduc- 
tion occurs on the subsystems unless they are irreducible. 
Irreducible imprimitive matrices can be written as 



/ C12 

c 23 












\ 





Ch-lh 

/ 



(66) 



where the blocks along the diagonal are square (second 
part of Perron Frobenius theorem) . The hth power of such 
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a matrix is block diagonal and the blocks are primitive 
|45| . so the stability analysis is again reduced. 

Physically, primitive matrices have the property that 
their powers are positive 5 (have no elements). From a 
physical perspective, primitive systems are thus "fully in- 
teracting" . This is in contrast to other nonnegative matri- 
ces which have zero blocks when raised to any power. 



C Further discussion of properties of A 

There is a correlation between an ill-conditioned A and 
a small eigenvalue gap. This is so because a matrix with 
a large k(X) is close to a matrix where A is repeated. In 
particular ;35 ; , there exists a matrix E such that A is a 
repeated eigenvalue of A + E and 

I A I 

IEI < 



However, k(A) may be small even if the gap is small. The 
relation between k(X) and the eigenvalue gap is shown in 
figure above. 

There is also a correlation between normality of A and 
a small k(A). When A is normal, that is, AA T = A T A, all 
of its eigenvectors are orthogonal and all the eigenvalues 
are perfectly conditioned. However, n(X) may be small in 
matrices which are far from normal. The relation between 
n(X) and the normality of A is shown in figure PHI 

Another way to characterize n is the relation 



i-£(u- 



^f)(v 



(67) 



where is the ith column of P (the right eigenvector 
corresponding to Xi) and ef is the «th row of P _1 (the 
left eigenvector corresponding to Ai) This relation is es- 
tablished by noting that e^e^efe^ = 1. It shows 

how v 2 is related to the angles between the eigenvectors. 
In particular, we see that for a normal matrix where the 
eigenvectors are orthogonal, v 2 = 1; but in general, the 
angular distribution of the eigenvalues is complicated. 

Finally, there is a correlation between | A2 1 — > and 
the variance a\ of the elements of A. Bounds for the sec- 
ond largest eigenvalue can be found in the case of row 
(or column) stochastic matrices, for example |45|: | A2 1 < 
min (l — X^i mm j Aij , maxj Aij — l) . This shows that, 
at least for stochastic matrices, a small variance u\ cor- 
responds to a large eigenvalue gap. This is shown to be 
true for all matrices in figure El Of course, the converse 
is not true; matrices with large a 2 A can also have a large 
eigenvalue gap, as also is shown in figure ITK1 



5 More exactly, the pth power of a nonnegative n x n prim- 
itive matrix A has no zero elements for all p > 7(A), where 
7(A) (the index of primitivity) is at most n 2 — 2n + 2, and 
usually much less|45)1. 



D Bounds on convergence of (|x| 2 ) 

In this section we apply the matrix 2-norm to determine 
two different bounds on the variance of the noise which, 
if satisfied, ensure the convergence of (|x| 2 ). These condi- 
tions are sufficient but by no means necessary. The second 
moment will of course never converge if the system does 
not converge in mean. We therefore take A < 1 in this 
section. 

The norm of a matrix is any function satisfying the 
regular properties of a vector norm and additionally the 
inequality |AB| < |A||B|. The matrix 2-norm correspond- 
ing to the usual Euclidean vector norm is 



|A| a = (p(AA*)) 



1/2 



(68) 



where p is the spectral radius. Note that for any norm, 
\A\ > |xA|/|x| = |A| for any eigenvalue A, so that in 
particular, 

A<|A|. (69) 

For ill-conditioned matrices, which includes those with ill- 
conditioned A, |A| is typically much larger than A 



D.l Bound on convergence of any moment 

We have |x* +1 |p < |A + B*| p |x*| p , so that 



<I*T> < 



t-i 



n<i A+BT 



.T=l 



where the expected value goes inside the product because 
the noise is white noise. (|x| p ) will thus converge for any 
p provided that (|A + B|) < 1 (we neglect the time super- 
script because the noise is stationary), or more usefully 



|B|)<1-|A| 



(70) 



using |A + B| < |A| + |B|. Since convergence of every 
moment is a much stronger condition than convergence 
of just the second moment, this bound is typically poor 
when applied to the second moment. 

We may estimate a lower limit for this bound for well- 
conditioned systems in the large n limit when the noise is 
UH (uncorrelated Bij all with the same variance b 2 ). We 
do so by using |A| > A (equation 169(1 and a result of |4"r>] 
that liminf |B| > 2b^/n almost surely in the large n limit, 
provided that the elements of B are mean i.i.d. and their 
moments do not grow too fast (which is satisfied for any 
reasonable noise). Thus 

Ai + 2b\fn < |A| + (|B|) 

and the condition ((70(1 on b for convergence at least weaker 
than 

t 1 -^ (71) 



b 2 



< 



An 



in the large n limit. That is to say, l(7U(l is more restrictive 
on b than 1(711) . For ill-conditioned systems, H71|l may not 
be accurate because |A| may be much larger than A. 
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Fig. 14. Scatter plots of |v| versus the Henrici number |AA — A A|, a measure of non-normality, for 10,000 randomly 
generated 5x5 matrices normalized so that A — 1. At left the elements were chosen from a normal distribution and only plotted 
if A was real. At right, the matrices are nonnegative primitive; a random number of elements were 0, and the nonzero elements 
were chosen from a uniform distribution. 
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Fig. 15. Scatter plots of eigenvalue gap A — | A2 1 versus standard deviation of the Aij for 10,000 randomly generated 5x5 
matrices, normalized so that A = 1. At left the elements were chosen from a normal distribution and only plotted if A is real. 
At right, the matrices are nonnegative primitive; a random number of elements were 0, and the nonzero elements were chosen 
from a uniform distribution. 



D.2 Second moment bound 

A different bound on the convergence of the second mo- 
ment in the case of UH noise can be found by applying 
the expected value before taking norms. We have 

(|x t | 2 ) = (|(A + B t )x*- 1 | 2 ) 

^(lAp+r^Klx'- 1 ! 2 ) 

where we have used the properties of the norm and the 
fact that the noise is UH, white and has mean 0. We thus 
have the convergence condition |A| 2 + nb 2 < 1, or 

62< (iz_W!2 

n 

for the convergence of (|a;| 2 ). Note that this condition is 
at least weaker than the condition 

b 2 K O^Af) (?2) 



because of l|69|l . Again, l|72|l may not be accurate for ill- 
conditioned A because |A| may be much larger than A. 
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